سناریوی مسئله:فرض کنید ماتریس کوواریانس سیستم به صورت زیر است که همبستگی بالای متغیرها (بهویژه بین همسایهها) را نشان میدهد. میانگین اولیه تمام متغیرها صفر است، اما مشاهده میکنیم که \(X_1 = 2\) رخ داده است.
تعریف ماتریس کوواریانس و مشاهداتتکهکد# ماتریس کوواریانس اصلی (Sigma)
Sigma <- matrix(c(1.0, 0.8, 0.6,
0.8, 1.0, 0.8,
0.6, 0.8, 1.0),
nrow=3, byrow=TRUE)
# میانگینهای اولیه
Mu <- c(0, 0, 0)
# پارامتر مشاهده شده (Observation)
x1_observed <- 2
print(Sigma)## [,1] [,2] [,3]
## [1,] 1.0 0.8 0.6
## [2,] 0.8 1.0 0.8
## [3,] 0.6 0.8 1.0محاسبات ماتریسی: بهروزرسانی باورهاقبل از شبیهسازی، باید توزیع شرطی \((X_2, X_3 | X_1=2)\) را محاسبه کنیم. این کار با استفاده از افراز ماتریس (Partitioning) انجام میشود.\[\mu_{new} = \mu_b + \Sigma_{21}\Sigma_{11}^{-1}(x_a - \mu_a)\]\[\Sigma_{new} = \Sigma_{22} - \Sigma_{21}\Sigma_{11}^{-1}\Sigma_{12}\]تکهکد# استخراج زیرماتریسها
Sigma_11 <- matrix(Sigma[1,1])
Sigma_12 <- Sigma[1, 2:3]
Sigma_21 <- matrix(Sigma[2:3, 1])
Sigma_22 <- Sigma[2:3, 2:3]
# محاسبه میانگینهای شرطی (بهروزرسانی شده)
# با توجه به اینکه Mu صفر است، فرمول ساده میشود
mu_new <- Sigma_21 %*% (1/Sigma_11) * (x1_observed - Mu[1])
# محاسبه کوواریانس شرطی (کاهش عدم قطعیت)
Sigma_new <- Sigma_22 - (Sigma_21 %*% (1/Sigma_11) %*% t(Sigma_12))
# نامگذاری برای خوانایی
rownames(mu_new) <- c("X2", "X3")
rownames(Sigma_new) <- c("X2", "X3")
colnames(Sigma_new) <- c("X2", "X3")
print("میانگینهای جدید (Updated Means):")## [1] "میانگینهای جدید (Updated Means):"## X2 X3
## [1,] 1.6 1.2## [1] "ماتریس کوواریانس جدید (Updated Covariance):"## X2 X3
## X2 0.36 0.32
## X3 0.32 0.64تفسیر نتایج محاسباتی:
- تغییر میانگین: مشاهده مقدار \(X_1=2\) باعث شد انتظار ما از \(X_2\) به ۱.۶ و از \(X_3\) به ۱.۲ افزایش یابد. این به دلیل همبستگی مثبت است.
- کاهش واریانس: واریانس \(X_2\) از ۱.۰ به ۰.۳۶ کاهش یافت. این یعنی مشاهده \(X_1\) اطلاعات زیادی درباره \(X_2\) به ما داده است و عدم قطعیت به شدت کم شده است.
شبیهسازی زنجیرهای (Chain Rule Simulation)حال که پارامترهای توزیع شرطی \((X_2, X_3)\) را داریم، از روش زنجیرهای برای تولید نمونه تصادفی استفاده میکنیم.این روش دقیقاً همان فرآیندی است که توابع کتابخانهای انجام میدهند، اما ما آن را به صورت دستی و گامبهگام (Manual Implementation) پیادهسازی میکنیم.استخراج پارامترهای دستی از ماتریس جدیدتکهکد# میانگینهای شرطی
1 واریانس و کوواریانس شرطی از ماتریس Sigma_new
اجرای الگوریتم شبیهسازیدر اینجا ابتدا \(X_2\) را تولید کرده و سپس \(X_3\) را بر اساس مقدار تولید شدهی \(X_2\) شبیهسازی میکنیم.تکهکدset.seed(123) # برای تکرارپذیری نتایج
2 گام ۱: شبیهسازی X2 از توزیع حاشیهای جدید
3 نکته: rnorm انحراف معیار میگیرد (جذر واریانس)
x2 <- rnorm(1, mean=mu_2, sd=sqrt(var_2))
# گام ۲: شبیهسازی X3 به شرط X2 (زنجیره دوم)
# محاسبه پارامترهای شرطی لحظهای برای X3
# ضریب رگرسیون X3 روی X2 در فضای جدید (beta)
beta <- cov_23 / var_2
# میانگین شرطی X3 با دانستن X2 تولید شده
cond_mean_3 <- mu_3 + beta * (x2 - mu_2)
# واریانس شرطی X3 (مقدار باقیمانده)
cond_var_3 <- var_3 * (1 - (cov_23^2 / (var_2 * var_3)))
# تولید نهایی X3
x3 <- rnorm(1, mean=cond_mean_3, sd=sqrt(cond_var_3))
# تجمیع نتایج در یک بردار
result_vector <- c(x1_observed, x2, x3)نتیجه یک بار اجرای شبیهسازی:
تکهکدprint(result_vector)همانطور که میبینید، مقدار تولید شده برای \(X_2\) (حدود ۱.۲۶) و \(X_3\) (حدود ۰.۸۸) حول میانگینهای محاسبه شده (۱.۶ و ۱.۲) نوسان دارند، اما دقیقاً برابر آنها نیستند که نشاندهنده ماهیت تصادفی فرآیند است.
تحلیل تصویری و اعتبارسنجیبرای اینکه رفتار سیستم را بهتر درک کنیم، شبیهسازی فوق را ۱۰۰۰ بار تکرار میکنیم تا توزیع احتمالاتی \(X_2\) و \(X_3\) را مشاهده کنیم.تکهکد# شبیهسازی انبوه برای تولید نمودار
n_sim <- 1000
sim_data <- data.frame(X2 = numeric(n_sim), X3 = numeric(n_sim))
for(i in 1:n_sim) {
# تکرار همان الگوریتم بالا
x2_loop <- rnorm(1, mean=mu_2, sd=sqrt(var_2))
beta_loop <- cov_23 / var_2
mean_3_loop <- mu_3 + beta_loop * (x2_loop - mu_2)
var_3_loop <- var_3 * (1 - (cov_23^2 / (var_2 * var_3)))
x3_loop <- rnorm(1, mean=mean_3_loop, sd=sqrt(var_3_loop))
sim_data$X2[i] <- x2_loop
sim_data$X3[i] <- x3_loop
}4 رسم نمودار توزیع توأم
ggplot(sim_data, aes(x=X2, y=X3)) +
# نقطه مرکزی تئوری
annotate("point", x=1.6, y=1.2, color="red", size=5, shape=18) +
geom_density_2d(color = "blue", size=1) +
geom_point(alpha=0.3, color="darkcyan") +
geom_vline(xintercept = 1.6, linetype="dashed", color="red") +
geom_hline(yintercept = 1.2, linetype="dashed", color="red") +
labs(
title = "فضای نمونهسازی شده برای (X2, X3) با شرط X1=2",
subtitle = "نقطه قرمز: میانگین تئوری محاسبه شده | خطوط آبی: منحنیهای چگالی",
x = "مقدار X2 (شبیهسازی شده)",
y = "مقدار X3 (شبیهسازی شده)"
) +
theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
تحلیل نمودار:ابر نقاط نشان میدهد که توزیع جدید همچنان بیضیگون (نرمال) است. مرکز ثقل دادهها دقیقاً روی نقطه \((1.6, 1.2)\) قرار گرفته است. همچنین کشیدگی بیضی در راستای قطر اصلی، نشاندهنده همبستگی مثبت باقیمانده بین \(X_2\) و \(X_3\) است.
خلاصه دستاوردهای شبیهسازی:
- قدرت روش شرطی:ما توانستیم بدون استفاده از توابع آمادهی جعبهسیاه (Black-box)، صرفاً با استفاده از قوانین جبر خطی و تولید اعداد تصادفی تکمتغیره، یک سیستم پیچیده چندمتغیره وابسته را شبیهسازی کنیم.
- انتشار اطلاعات (Information Propagation): مشاهده کردیم که چگونه دیدن مقدار \(X_1\) نه تنها میانگین همسایه نزدیکش (\(X_2\)) را تغییر داد، بلکه اثر آن تا متغیر دورتر (\(X_3\)) نیز منتشر شد، اگرچه شدت اثر (همبستگی) کمتر بود.
-
کدنویسی منعطف: الگوریتم پیادهسازی شده (محاسبه
betaو پارامترهای لحظهای) پایهای برای روشهای پیشرفتهتر مانند Gibbs Sampling در آمار بیزین است.
کاربرد عملی:
این روش در مهندسی مالی (برای پیشبینی قیمت داراییهای وابسته)، هواشناسی (پیشبینی وضعیت نقاط جغرافیایی مجاور) و پردازش سیگنال (فیلتر کالمن) کاربرد حیاتی دارد.